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We present a global hybrid meta-generalized gradient approximation (meta-GGA) with three 
empirical parameters, as well as its underlying semilocal meta-GGA and a meta-GGA with only 
one empirical parameter. All of them are based on the new meta-GGA resulting from the under- 
standing of kinetic-energy-density dependence [J. Chem. Phys. 137, 051101 (2012)]. The obtained 
functionals show robust performances on the considered molecular systems for the properties of 
heats of formation, barrier heights, and noncovalent interactions. The pair-wise additive dispersion 
corrections to the functionals are also presented. 



I. INTRODUCTION 

The Kohn-Sham (KS) density functional theory (DFT) 
[IH3] is one of the most widely used electronic struc- 
ture theories for atoms, molecules and solids in vari- 
ous areas of physics, chemistry and molecular biology 
due to its computational efficiency and useful accu- 
racy. It simplifies a many-electron wave-function prob- 
lem to an auxiliary one-electron problem, with only its 
exchange-correlation part carrying the many-electron ef- 
fects to be approximated in practice. Among numer- 
ous exchange-correlation approximations, the local spin 
density approximation (LSDA) [TJ H] [5], the standard 
Perdew-Burke-Ernzerhof (PBE) generalized gradient ap- 
proximation (GGA) [6 , and the Becke-3-Lee- Yang-Parr 
(B3LYP) [7HTU] hybrid GGA dominate the user market 
of DFT QT] . The former two are efficient semilocal func- 
tionals widely used for extended systems, while the latter 
is a computationally more expensive nonlocal functional 
that hybridizes a GGA with the exact exchange energy 
and is popular for finite systems. At the semilocal level, 
however, meta-GGA (MGGA) is the highest rung of the 
so-called Jacob's ladder in DFT [12] and potentially the 
most accurate one [1.3!, which can also serve as a better 
base for hybridizing with the exact exchange energy. 

Semilocal approximations (e.g., Refs. [4H6] H4HT6] ) of 
the form 

-SxcK* n i] = J d 3 rne S x C (n t , 714., Vra t , Vr^, r t , tj.) (1) 

require only a single integral over real space and so are 
practical even for large molecules or unit cells. In Eq. 
([T]), n-j- and are the electron densities of spin ^ and 
4,, respectively, Vn^j the local gradients of the spin den- 
sities, t^j. = ^2 k IVV'fct,^ 2 /2 the kinetic energy densi- 
ties of the occupied KS orbitals yj krT of spin a {a =t,l), 
and e^c the approximate exchange-correlation energy per 
electron. All equations are in atomic units. In addi- 
tion to n-f and nj,, the only ingredients in LSDA, GGAs 



also use the density gradients; MGGAs additionally in- 
clude the kinetic energy densities r a . With the en- 
coded information of shell structures from the kinetic en- 
ergy densities, MGGAs can distinguish different orbital- 
overlap regions |17j . and thus deliver simultaneous accu- 
rate ground-state properties for molecules, surfaces, and 
solids - which couldn't be obtained with a GGA or LSDA 
[13] . Computationally, MGGAs are not much more ex- 
pensive than LSDA or GGA [T51 [T^] . In computations 
for molecules containing transition-metal atoms, the Tao- 
Perdew-Staroverov-Scuseria (TPSS) iT5] MGGA is only 
30% slower [19] than PBE. 

To understand the performance of MGGAs, a re- 
cent study [17] investigated the effect of r-dependence 
on MGGAs through the inhomogeneity parameter a = 
(r - T W )/r UDii . Here, r = 52 a r a , t w = ||Vn| 2 /n is 
the von Weizsacker kinetic energy density, and r unlf = 

(37r 2 ) 2 / 3 n 5 / 3 is the kinetic energy density of the uni- 
form electron gas (UEG). t w is a lower bound on r with 
t w — t only for a single-orbital region [20] . Therefore, 
a ^ measures locally how far n(r) deviates from being 
of single-orbital character on the scale of the UEG, and 
thus characterizes the extent of orbital overlap. Besides 
the ability of a to distinguish different orbital-overlap 
regions, Ref. [T7| further showed that the effect of the 
a-dependence on MGGAs is qualitatively equivalent to 
that of the dependence on the reduced density gradient 
s = I Vn|/[2(37r 2 ) 1 / 3 n 4 / 3 ], another dimensionless inhomo- 
geneity parameter that measures how fast and how much 
the density varies on the scale of the local Fermi wave- 
length 27T jkp with kp = (37T 2 n) 1 / 3 . The s-dependence is 
well understood at the GGA level [SHU], which MGGAs 
inherit [HHIZ]- 

As a result of the study on the r-dependence [17], a 
simple exchange functional, where the a-dependence is 
disentangled from the s-dependence, was constructed as 
an interpolation between the single-orbital regime (a — 
0) and the slowly- varying density regime (a « 1) [17j . 
which we will discuss more in Section |III[ In analogy to 
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the monotonically increasing s— dependence that GGAs 
are often designed to have to favor more inhomogeneous 
electron densities, the simple exchange functional has a 
monotonically decreasing a— dependence to favor more 
the single-orbital regions. When combined with a variant 
of the PBE correlation [16j OH, the resulting MGGAs 
(denoted as MGGA_MS with MS standing for "made 
simple") |17j performs equally well for atoms, molecules, 
surfaces, and solids, with an overall performance that is 
comparable to the sophisticated revised TPSS (revTPSS) 
MGGA [TB]. Moreover, MGGAJVIS yields excellent 
binding energies for the W6 water clusters [T7], even 
better than those from the highly parametrized M06-L 
MGGA [21] whose training sets include noncovalent in- 
teractions (hydrogen bonds and van der Waals interac- 
tions). Further tests on general main group thermochem- 
istry, kinetic, and noncovalent interactions [2 2) showed 
that MGGAJVIS, though not as good as M06-L, system- 
atically improves the descriptions for the noncovalent in- 
teractions over conventional semilocal ftmctionals, e.g., 
PBE, TPSS, and revTPSS. 

By including training sets of noncovalent interactions, 
the M06-L MGGA was trained to capture medium-range 
exchange and correlation energies that dominate equilib- 
rium structures of noncovalent complexes |21j . For exam- 
ple, M06-L delivers good binding energies for the S22 set 
that is dominated by noncovalent interactions and will be 
discussed more in Sec. |IV[ good structural and energetical 
properties of layered materials bound by van der Waals 
interactions (e.g., graphite, hexagonal boron nitride, and 
molybdenum disulfide |23j). and good adsorption ener- 
gies of aromatic molecules adsorbed on coinic metal sur- 
faces [23]. M06-L represents one of the developments of 
MGGAs that are fitted to a variety of training sets with 
numerous parameters to try to improve different proper- 
ties together [2JJE5, ^6 ■ However, it should be stressed 
that too many fitting parameters could cause problems, 
e.g., convergence of calculations and the related smooth- 
ness of binding curves |27j , and the prediction for systems 
far from training sets might not be realistic |28j . 

Semilocal (si) approximations can be reasonably ac- 
curate for the near-equilibrium and compressed ground- 
state properties of ordinary matter, where neither strong 
correlation nor long-range van der Waals interaction is 
important. However, stretched bonds that arise in tran- 
sition states of chemical reactions or in the dissociation 
limits of some radical or heteronuclear molecules, etc., 
require full nonlocality |29| . The global hybrid ftmction- 
als discussed later can provide such needed nonlocality. 
Long-range van der Waals interactions, originated from 
spontaneous correlations between two distant electron 
densities, also require full nonlocality |30) . which is dif- 
ferent from the previous one [31) . To obtain long-range 
van der Waals interactions, there are various DFT-based 
dispersion techniques, including the post-DFT empirical 
pairwise potential corrections [32TI39] and the van der 
Waals density functionals (vdW-DFs) [Ml gD] 2JJ , which 
are based on the pairwise additivity and problematic for 



metallic systems :42, l.'i . See Ref. HH for an overview of 
the progess on DFT-based dispersion techniques. Con- 
ventional semilocal functionals, e.g., the widely used PBE 
GGA, even lack the ability to describe noncovalent inter- 
actions near equilibrium, while the M06-L MGGA and 
MGGAJVIS as well as its variant proposed here improve 
over PBE. 

The global hybrid (gh) idea, due to Becke [7] 25] , intro- 
duces some full nonlocality into the calculation but only 
at the level of £"° xact , which can be evaluated semiana- 
lytically from the Kohn-Sham orbitals in some computer 
codes. In its simplest (one parameter) version [46) : 

Eg = aEl™« + (l-a)E% + Ef (2) 

where the exact-exchange mixing parameter a takes an 
empirical value in the range < a < 1 . A rough but con- 
vincing argument |47j for why Eq. ^ works is: semilocal 
exchange-correlation typically overestimates atomization 
energies and underestimates energy barriers of chemi- 
cal reactions, while exact exchange without correlation 
makes errors of the opposite sign, so that mixing of the 
two will yield better atomization energies and barriers 
than either alone. Therefore, values of a vary with the 
choice of semilocal functionals when fitting Eq. ([2| to 
formation enthalpies. The more a semilocal functional 
overbinds molecules, the larger the value of a is. Csonka 
et al. E7] found a to be 0.60 for the PBEsol GGA, 0.32 
for the PBE GGA, and 0.1 for the TPSS and revTPSS 
MGGAs. Improving the underlying semilocal functional 
reduces the value of a needed to fit the atomization ener- 
gies or enthalpies of formation, which however can worsen 
the barrier heights. Some global hybrids [3H1 EH] thus 
employ highly fitted semilocal parts intended not to be 
accurate by themselves but to work well with a fraction 
of exact exchange. 

The popular B3LYP JTHTD] hybrid GGA has a some- 
what more complicated form and reads: 

£B3LYP = a £exa.t + (1 _ a) £LSDA + 6A £B88 + 

+ (1 - c)£ f VWN3 + c£ c LYP (3) 

where E^ WN3 [5U] is the LSDA correlation energy in 
the random phase approximation which does not yield 
the correct uniform electron gas limit. In the equa- 
tion, a=0.2, b=0.72, and c=0.81. If we use the simi- 
lar idea as in Eq. ^ to define the underlying semilo- 
cal functional for B3LYP by replacing E° xact with £| 88 , 
we end up with a GGA, denoted as BLYP*, E^ YP = 

£LSDA + (fl + b)AE B8S + (j_ _ C )£VWN3 + ^LYP BLYP* 

differs from BLYP in reducing the gradient correction of 
the B88 exchange to (a+b)=92% and mixing the VWN3 
and LYP correlations. This difference confirms the rough 
argument of the previous paragraph. Typically, BLYP 
underbinds molecules, and inclusion of exact exchange 
worsens the underbinding. The reduction of the gradi- 
ent correction of the B88 exchange and the mixture of 
the VWN3 and LYP correlations then turn BLYP* into 
an overbinding functional suitable for hybridization with 
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the exact exchange. The idea of tuning the gradient cor- 
rection of B3LYP can be transfered to the construction 
of our hybrid MGGA, where both the s— dependence and 
the a— dependence are tuned. 

In the next section, we will present the computational 
details, followed by the constructions of semilocal and 
hybrid MGGAs in Section |III| Results and discussions 
will be given in Section |IV| And we give the conclusions 
in Section IVl 



II. COMPUTATIONAL DETAILS 

The functionals described here were implemented in 
the development version of the Gaussian electronic 
structure program [5TJ . All calculations employ the fully 
uncontracted 6-311+- {-G(3d/, 3pd) basis set jSHES] to ob- 
tain benchmark quality results. We used the UltraFine 
grid with 99 radial shells and 590 angular points for 
numerical integration of the DFT XC potential. The 
Gaussian program was used for all calculations pre- 
sented in this study. 

Different training sets were used to adjust empirical 
parameters. For the AE6 and BH6 tests sets [53] we 
employed reference data from highly accurate quantum- 
chemical calculations |55j . Experimental reference val- 
ues were used for G2/97 (148 molecules) jSH [ST] and 
BH42/03 (21 forward and reverse hydrogen transfer bar- 
rier heights) [55]. The fitted parameters would not 
change much if the highly accurate quantum-chemical 
reference data from Ref. [55] would have been used. 

The superset of G2/97 and G3-3 (75 molecules) [5U] 
(a total of 223 molecules comprising the G3/99 test 
set) [ST], the S22 test set for weak interactions [52], as 
well as the barrier height test sets HTBH38/04 (19 for- 
ward and reverse hydrogen transfer barrier heights) [58] 
and NHTBH38/04 (19 forward and reverse non-hydrogen 
transfer barrier heights) [63] are employed for assessment 
of the fitted functionals. The employed geometries and 
reference values are available from the Supporting Infor- 
mation of Refs. [B2] and The counterpoise correction 
[65] was used to reduce the basis set superposition er- 
rors for calculations of weak interactions. A vibrational 
scale factor of 0.9854 was used for calculations of heats 
of formation on the B3LYP/6-31G(2d/, p) level of theory 
as recommended in Ref. 1661 

The performance of our fitted functionals will be com- 
pared to related and popular density functionals: PBE 
0, TPSS [IS], revTPSS [JJ], M06-L [UJ, MGGAJV1S0 
QI], PBEh 06], TPSSh [67], revTPSSh [16], B3LYP 0- 
HD], BLYP [TO], and M06 [68]. Calculations of open- 
shell species were carried out via spin-unrestricted for- 
malisms. Errors are reported as calculated — reference. 



III. CONSTRUCTION OF SEMILOCAL AND 
HYBRID META-GGAS 

The semilocal exchange energy of a spin-unpolarized 
density can be written as: 

E*[n] = J d s rne x ^{n)F x {p,a). (4) 

Here, e" nlf (n) = — j^i^jr) 1 ^ l r s is the exchange energy 
per particle of a UEG with r s — (47m/3) -1 / 3 , p = s 2 , 
and F x the enhancement factor with F x = 1 for LSDA. 
The expression for the spin-polarized case follows from 
the exact spin scaling [55]. In Ref. [TTJ a simple ex- 
change enhancement factor that disentangles the a- and 
p-dependence was proposed, 

F^ (p, a) = Fi (p) + /(at) [F° (p) - F^ (p)] , (5) 

where F*(p) = F x nt (p, a = l) = l + «; - k/(1 + pP^p/n) 
andF°(p) = F x nt (p,a = 0) = l+n-n/(l + (p GE p+c)/n). 
F x nt (p,a) interpolates between F x (p) and F x (p) through 
f(a). Here, in order to tune the a dependence, we 
add a parameter 6 in the denominator of the original 
/(a), which reads now /(a) = (1 - a 2 ) 3 /(l + a 3 + 6a 6 ). 
For a slowly varying density where awl, /(a) is con- 
trolled by its numerator and we recover the second or- 
der gradient expansion with the first principle coefficient 
^ge = 10/81 [70] as in Ref. 17, The parameter c is de- 
termined for each k to reproduce the exchange energy of 
the hydrogen atom, where a — 0. Since 6a 6 is negligible 
for < a < 1 and the numerator of f(a) modulates the 
behavior of F x nt (p, a) for a rj 1,6 mainly controls the be- 
havior for a > 1. On the other hand, k mainly controls 
the behavior of F x nt (p, a) for large s. In the original 
form [T7], 6 = 1, and k — 0.29 was determined by the 
Hartree-Fock exchange energy of the 12-noninteracting- 
electron hydrogenic anion with nuclear charge Z=l |71j . 
Fig. 2 of Ref. 17! showed that the shell regions of the 
hydrogenic anion are associated with a < 1 and the in- 
tcrshcll regions with a > 1, which thus can provide in- 
formation to guide the exchange functional from a ~ 1 
to a — > oo. Here, we release the constraint of the hydro- 
genic anion and rely on some molecular training sets to 
provide information to determine both 6 and n by fitting 
them to heats of formation and barrier heights. For the 
correlation part, we use the variant of the PBE correla- 
tion PHOT]. 

When we plug the above exchange and correlation 
parts into Eq. ([2]), we obtain a global hybrid MGGA 
with three parameters, k, 6, and a, which control the de- 
pendence of the global hybrid on the density gradient, the 
kinetic energy density, and the exact exchange, respec- 
tively. In the next section, we describe the parametriza- 
tion procedure and the performance of our best function- 
als. 



4 




0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 



FIG. 1: The AMAE surface for the semilocal MGGA_MSi 
variants with errors in kcal/mol shows a broad range of rather 
small errors 



IV. RESULTS AND DISCUSSIONS 

To get a broad overview about the errors which the 
semilocal functional class MGGA_MS makes using dif- 
ferent values of n and b, we define the average mean 
absolute error (AMAE) using the mean absolute error 
(MAE) values from the AE6 and BH6 test sets: 

AMAE = 0.5 [MAE(AE6) + MAE(BH6)] , (6) 

and present the plotted surface of AMAE values in figure 
[T] The figure shows that there are two useful possibili- 
ties of fitting a semilocal functional of our MGGA_MS 
form: First, we can design a one-parameter functional 
by setting the parameter b to one and fitting n to min- 
imize a chosen error measure. Second, both parameters 
can be fitted independently. The second choice will yield 
a lower AMAE value. For the independent fit, we also 
include the exact exchange admixture parameter a into 
the fitting procedure. The corresponding semilocal func- 
tional will be obtained by setting a to zero while keeping 
the other global hybrid parameters fixed. 

Our chosen error measure contains the mean absolute 
errors (MAE's) of the heats of formation of the G2/97 
data set and barrier heights of the BH42/03 test sets: 

SW = 0.5 [MAE(G2/97) + MAE(BH42/03)] . (7) 

The resulting parameter sets, ME, and MAE values are 
shown in table [T] and compared to other popular semilo- 
cal and hybrid functionals. MGGAJV1S0 is the original 
MGGA_MS functional from Ref. [17] with no parameter 
fitting to training sets while MGGA_MS1 is our semilo- 
cal one-parameter functional. MGGA_MS2h is the new 
hybrid functional and MGGAJMS2 its semilocal comple- 
ment. Our hybrid MGGA_MS2h has a rather low ex- 
act exchange admixture of 9% as expected when a good 
semilocal functional is hybridized [47] . Comparing the 



parameter set of MGGA_MS2 with the AMAE surface 
in figure [T] shows that its parameter set is quite close 
to the global minimum of the AMAE, although it has 
not been refitted after exact exchange admixture was re- 
duced from 9% to 0%. The functional BLYP* is not com- 
monly used, but it is the underlying semilocal functional 
of B3LYP. ft is included in our tables for comparison pur- 
poses only. As shown in tables [T] and \H\ BLYP* turns the 
underestimations of BLYP for atomization energies into 
overestimations. From all semilocal functionals in ta- 
ble JTJ MGGAJMS1 performs best for G2/97 while M06-L 
performs best for BH42/03. MGGAJMS0, MGGAJVISl, 
and MGGAJV1S2 functionals are the next best perform- 
ers for the BH42/03 test set, which behave similarly 
despite very different parameter sets. For the hybrids, 
B3LYP performs best for G2/97 and M06 for BH42/03. 
Our global hybrid MGGA_MS2h performs similarly to 
B3LYP and PBEh for the BH42/03 test set and similar 
to PBEh for the G2/97 test set. Note that PBE, TPSS, 
rcvTPSS, and MGGA_MS0 are nonempirical functionals, 
while the number of empirical parameters - as indicated 
in the parentheses - of MGGAJVISl (1), MGGA.MS2 
(2), TPSSh (1), revTPSSh (I), and MGGA_MS2h (3) 
are an order of magnitude smaller than those of M06-L 
and M06. 

MGGAJVISl performs better than MGGAJVIS2 and 
MGGA_MS0 for all molecular test sets discussed so far. 
This and figure [l] indicate that one can construct another 
MGGA_MS functional which performs even better than 
MGGA_MSI. However, figure [T] suggests that the im- 
provement over MGGA_MS1 will not be significant. We 
found that k = 0.514, c = 0.14352, and b = 2.0 is close to 
optimal for the error measure SW, but the minimum in 
parameter space is rather shallow and broad. With MAE 
values for G2/97 of 4.5 kcal/mol and for BH42/03 of 5.6 
kcal/mol, this variant would improve over MGGAJVISl 
only slightly for G2/97 and not at all for BH42/03 on the 
expense of an additional empirical parameter. Thus, we 
drop this variant from now on. 



To test the deviation from the exact exchange energy 
of the 12-noninteracting-electron hydrogenic anion, wc 
calculate the exchange energies for the new parameter 
sets of (b, k), which are -1.8562 Ha for (b=l, k = 0.404) 
of MGGA_MS1 and -1.8558 Ha for (b=4, k = 0.504) 
of MGGA_MS2, respectively. In comparison with the 
exact one, -1.8596 Ha, the deviations are small, which 
indicates that the information of the hydrogenic anion 
is preserved by Eq. [5] and complemented by that from 
training sets. Figure [2] shows the exchange enhancement 
factors for different a and s values of our new functionals 
and compares them to the revTPSS functional. Figure 
shows revTPSS has the order of limits anomaly at 
1 s and small a [151 116] , which can be eliminated by 
regularizing revTPSS at these regions |72) . The family of 



2(a) 



5 



TABLE I: Selected parameters, mean errors (ME), and mean absolute errors (MAE) of semilocal and hybrid functionals for 
the G2/97 and BH42/03 test sets. The ME and MAE values are in kcal/mol. 





a 


K 


c 


b 


G2/97 


BH42/03 












ME 


MAE 


ME MAE 


PBE 


0.00 


0.804 




— 


-16.1 


16.9 


-9.7 


9.7 


BLYP 


0.00 






- 


-0.6 


7.4 


-8.0 


8.0 


BLYP* 


0.00 






— 


-8.1 


9.0 


-8.4 


8.4 


TPSS 


0.00 


0.804 




- 


-5.7 


6.4 


-8.4 


8.4 


revTPSS 


0.00 


0.804 




- 


-4.7 


5.7 


-7.7 


7.7 


M06-L 


0.00 








-3.0 


5.1 


-4.4 


4.4 


MGGAJVISO 


0.00 


0.29 


0.28771 


1.0 


-0.6 


6.8 


-5.9 


6.0 


MGGA.MS1 


0.00 


0.404 0.18150 


1.0 


2.3 


4.9 


-5.5 


5.6 


MGGAJV1S2 


0.00 


0.504 0.14601 


4.0 


-0.8 


5.1 


-5.8 


5.8 


PBEh 


0.25 


0.804 






-2.6 


5.0 


-4.7 


4.7 


TPSSh 


0.10 


0.804 






-1.9 


4.4 


-6.6 


6.6 


revTPSSh 


0.10 


0.804 






-1.2 


4.5 


-6.0 


6.0 


B3LYP 


0.20 








0.9 


3.1 


-4.7 


4.7 


M06 


0.27 








-3.4 


4.1 


-2.3 


2.3 


MGGA_MS2h 0.09 0.504 0.14601 4.0 


2.7 


4.9 


-4.5 


4.6 



TABLE II: Mean errors (ME), and mean absolute errors (MAE) of semilocal and hybrid functionals for the G3-3, G3/99, 
HTBH38, and NHTBH38 test sets in kcal/mol 



G3-3 G3/99 HTBH38/04 NHTBH38/04 





ME 


MAE 


ME 


MAE 


ME 


MAE 


ME 


MAE 


PBE 


-32.6 


32.6 


-21.6 


22.2 


-9.7 


9.7 


-8.5 


8.6 


BLYP 


12.5 


14.1 


3.8 


9.7 


-7.9 


7.9 


-8.7 


8.8 


BLYP* 


-5.7 


8.6 


-7.3 


8.8 


-8.4 


8.4 


-8.7 


8.8 


TPSS 


-6.4 


6.5 


-6.0 


6.5 


-8.2 


8.2 


-9.0 


9.1 


revTPSS 


-4.2 


4.7 


-4.5 


5.4 


-7.4 


7.4 


-9.3 


9.3 


M06-L 


-4.2 


7.4 


-3.4 


5.8 


-4.5 


4.5 


-3.2 


3.7 


MGGAJVIS0 


-5.9 


11.8 


-2.4 


8.5 


-5.6 


5.7 


-6.4 


6.9 


MGGA_MS1 


3.2 


5.2 


2.6 


5.0 


-5.4 


5.5 


-6.3 


6.7 


MGGAJV1S2 


-4.1 


7.6 


-1.9 


6.0 


-5.7 


5.7 


-7.0 


7.3 


PBEh 


-9.6 


10.4 


-4.9 


6.8 


-4.7 


4.7 


-3.2 


3.7 


TPSSh 


-1.0 


3.5 


-1.6 


4.1 


-6.4 


6.4 


-6.9 


7.0 


revTPSSh 


0.6 


3.6 


-0.6 


4.2 


-5.8 


5.8 


-7.1 


7.2 


B3LYP 


7.9 


8.2 


3.3 


4.8 


-4.5 


4.6 


-4.6 


4.7 


M06 


-4.4 


5.4 


-3.7 


4.5 


-2.3 


2.3 


-1.7 


2.2 


MGGA_MS2h 


1.5 


5.1 


2.3 


4.9 


-4.4 


4.4 


-5.3 


5.7 



MGGAJMS doesn't have such a problem by construction. 
The enhancement factors of different MGGAs considered 
here are on top of each other around aw 1 at the curves 
of s = 0.001, due to the constraint of the second order 



gradient expansion. Figure [2(b) shows that each MGGA 
is insensitive to a for large s. For the a = curves, which 
is modulated by the constraint of the exchange energy of 
the hydrogen atom, F x gets smaller for small s regions 
when a larger k is used. The behavior of these MGGAs at 
large s is determined by the value of k. The larger k, the 
larger enhancement factor at large s. Due to the qualita- 
tive equivalence between the s— and a— dependence [17] . 
we observe for different MGGAs that the faster the in- 
crease of F x with s, the slower the decrease of F^ with 
a. 

It is not too surprising that our functionals fitted to 
G2/97 and BH42/03 are among the best performers for 



these test sets. Table [IT] shows the performance for the 
heats of formation of the test sets G3-3 and G3/99 as well 
as the barrier heights of the data sets HTBH38/04 and 
NHTBH38/04. The Minnesota functionals perform best 
for the barrier heights, which is not too surprising as both 
barrier height test sets are contained in the training set of 
M06 and M06-L. The next best semilocal functional for 
both barrier height test sets is MGGA_MS1. Among the 
global hybrids, PBEh and B3LYP perform better than 
MGGA_MS2h for the NHTBH38/04 test set but not for 
the HTBH38/04 test set which was part of the training 
set. TPSSh performs best for G3-3 and G3/99 within our 
choice of global hybrids. Interestingly, all global hybrids 
in tab le[n] (except PBEh) perform very similar on average 
for the G3/99 test set. Among the semilocal functionals 
revTPSS performs best for G3-3 and MGGA.MSl for 
G3/99. 
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Aside from heats of formation and barrier heights, we 
would like to test our functionals for a property which 
was not part of the training set. We choose the weak in- 
teractions of the S22 test set. Table [TTT| presents the ME 
and MAE values for the S22 test set and its subsets WI8 
(including 8 dispersion-bound complexes), MI7 (includ- 
ing 7 mixed complexes), and HB7 (including 7 hydrogen- 
bonded complexes). Before discussing the results, we 
would like to test the convergence of MGGAs with re- 
spect to the integration grid for dispersion-bound com- 
plexes. Johnson et al [73) found that very large integra- 
tion grids (much finer than the UltraFine) are required 
to remove oscillations in potential energy surfaces (PES) 
for dispersion-bound complexes for most of the MGGAs 
they studied. For the chosen simple dispersion-bound 
NeAr, we confirmed that the M06L MGGA, which was 
in Johnson's test list, needs a very large grid to yield a 
smooth PES. However, our MGGAs converge much faster 
for the total energies than M06L does, and yield con- 
verged smooth curves with the UltraFine grid. The Min- 
nesota functionals, which have been fitted for weak inter- 
actions, perform best for all these data sets of weak in- 
teractions. Among semilocal functionals, MGGA_MS0 is 
the second best performer for the S22 set and its subsets. 
MGGAJV1S2 on average is comparable to MGGAJV1S0 
for the WI8 and MI7 subsets, but worse for the HB7 sub- 
set, while MGGA_MS1 loses the accuracy for the weak 
interactions compared to the other two in the MGGA_MS 
family. As expected for WI8 and MI7, where dispersion 
interactions are important, PBE, TPSS, and revTPSS 
yield large MAEs, which are worse than the other semilo- 
cal functionals considered here. 



Ref. [44] presents a collection of different techniques 
for functionals sorted by the way weak interactions are 
accounted for: step 1 (DFT-D 37] , etc), step 2 (DFT- 
D3 [35], vdW(TS) [39], etc), step 3 (long-range density 
functionals [HOI EB HI]), and step 4 and higher (RPA [74], 
etc) . Representative functionals of each group are tested 
for the S22 test set in Ref. HH Note that the basis sets 
differ compared to our calculations, but the results can be 
compared roughly. MGGAJV1S0 and MGGA_MS2 yield 
results that are close to the vdW-DF [3D] which is on 
step 3 of this classification, while MGGA_MS1 fits best 
to the ground (below step 1). It should be noted that the 
classification is methodological and not based on perfor- 
mance as some functionals in step 1 and 2 outperform 
some funtionals on step 3 in Ref. HU Table |IV| gives 
results of the MGGAJV1S functionals with the D3 cor- 
rection of Grimme et al [38 . The parameters in the D3 
correction term were optimized by fitting to the S22 data 
set for each functional. We found that the optimization 
is insensitive with respect to the parameter s&, and thus 
set sg = 0. Then, there is only one fitting parameter s r fi 
left in the D3 correction term. s r 6 , which is the scal- 
ing factor of cutoff radii and determines the range of the 




0.5 1 1.5 2 2.5 3 3.5 



s 

(b) 

FIG. 2: (a) Exchange enhancement factors vs. a for different 
s. The solid curves are for s = 0.001 and dashed for s = 1.001. 
(b) Exchange enhancement factors vs. s for different a. In 
addition to the PBE and PBEsol curves, the solid curves are 
for a — and dashed for a = 1. 



dispersion correction inversely [35], is 1.15, 1.05, 1.14, 
and 1.14 for MGGA_MS0, MGGA_MS1, MGGAJVIS2, 
and MGGA_MS2h, respectively. All the four MGGAJVIS 
functionals with the D3 corrections deliver for the S22 
data set similar MAEs around 0.3 kcal/mol, which are 
among the best of the D3-corrected functionals tested 
in Ref. [351 The smaller s rfi for MGGAJVISl indicates 
that MGGA_MS1 captures less intermediate-range weak 
interactions than the other three do, as also shown in 
Table IIII1 Table IIVI also shows that the D3 corrections 
have noticeable effects on the formation energies and 
the barrier heights. MGGAJV1S1+D3 is the best among 
the three MGGAs, while the hybrid MGGAJVlS2h+D3 
improves further the formation energies and the barrier 
heights. 

The functional pairs of global hybrid and its corre- 
sponding semilocal complement perform very similarly 
for weak interactions on average. The pair of BLYP* 
and B3LYP is the outlier, which might be related to the 
fact that B88 is too repulsive for rare gas dimers even 
when comparing to the exact exchange |75j . This gen- 
eral trend suggests that the performance of a hybrid of 
a semilocal functional inherits its performance for weak 
interactions from its corresponding pure semilocal func- 
tional. A possible explanation is that weak interactions 
are mainly correlation effects where exact exchange does 
not contribute. On the other hand, heats of formation, 
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TABLE III: Mean errors (ME), and mean absolute errors (MAE) of semilocal and hybrid functionals for the S22 test set and 
its subsets in kcal/mol 



WI8 


MI7 


HB7 


S22 


ME MAE 


ME MAE 


ME MAE 


ME MAE 



PBE 


4.8 


4.8 


2.0 


2.0 


1.1 


1.1 


2.8 


2.8 


BLYP 


7.7 


7.7 


3.8 


3.8 


3.2 


3.2 


5.0 


5.0 


BLYP* 


7.3 


7.3 


3.6 


3.6 


2.8 


2.8 


4.7 


4.7 


TPSS 


6.0 


6.0 


2.7 


2.7 


1.9 


1.9 


3.7 


3.7 


revTPSS 


4.3 


4.3 


2.4 


2.4 


2.0 


2.0 


2.9 


2.9 


M06-L 


1.0 


1.0 


1.0 


1.0 


0.7 


0.7 


0.9 


0.9 


MGGAJVISO 


3.0 


3.0 


1.2 


1.2 


0.8 


0.8 


1.7 


1.7 


MGGAJV1S1 


3.9 


3.9 


1.7 


1.7 


1.8 


1.8 


2.5 


2.5 


MGGAJVIS2 


3.0 


3.0 


1.3 


1.3 


1.4 


1.4 


1.9 


1.9 


PBEh 


4.6 


4.6 


1.8 


1.8 


0.8 


0.9 


2.5 


2.5 


TPSSh 


5.8 


5.8 


2.6 


2.6 


1.8 


1.8 


3.5 


3.5 


revTPSSh 


5.0 


5.0 


2.3 


2.3 


1.8 


1.8 


3.1 


3.1 


B3LYP 


6.5 


6.5 


3.0 


3.0 


2.0 


2.0 


4.0 


4.0 


M06 


1.4 


1.4 


0.9 


0.9 


0.7 


0.7 


1.0 


1.0 


MGGA_MS2h 


3.0 


3.0 


1.3 


1.3 


1.2 


1.2 


1.9 


1.9 



TABLE IV: Mean errors (ME), and mean absolute errors (MAE) of semilocal and hybrid functionals with the D3 correction |38| 
for the S22 test set and its subsets in kcal/mol. All functionals use s$ — and sg = 1. See Ref. 1381 for the definitions of se,sg, 
and s r ,6- 

MGGA.MS0+D3 MGGA.MS1+D3 MGGA.MS2+D3 MGGA_MS2h+D3 





Sr,6 


= 1.15 




= 1.05 


<Sr,6 " 


= 1.14 


$r,6 


= 1.14 


test set 


ME 


MAE 


ME 


MAE 


ME 


MAE 


ME 


MAE 


G2/97 


-0.9 


7.0 


1.8 


5.0 


-1.1 


5.3 


2.3 


4.9 


BH42/03 


-6.1 


6.1 


-5.9 


5.9 


-6.0 


6.0 


-4.7 


4.7 


G3-3 


-7.4 


13.1 


0.6 


5.7 


-5.6 


9.0 


-0.0 


5.8 


G3/99 


-3.1 


9.1 


1.4 


5.3 


-2.6 


6.6 


1.6 


5.2 


HTBH38/04 


-5.8 


5.9 


-5.8 


5.8 


-6.0 


6.0 


-4.6 


4.6 


NHTBH38/04 


-6.5 


7.0 


-6.5 


6.9 


-7.1 


7.4 


-5.4 


5.8 


WI8 


-0.0 


0.4 


-0.2 


0.4 


-0.2 


0.3 


-0.1 


0.3 


MI7 


-0.3 


0.3 


-0.3 


0.3 


-0.2 


0.3 


-0.2 


0.3 


HB7 


-0.3 


0.4 


0.2 


0.2 


0.2 


0.3 


0.1 


0.2 


S22 


-0.2 


0.3 


-0.1 


0.3 


-0.0 


0.3 


-0.1 


0.2 



atomization energies, and barrier heights are significantly 
affected by correlation and exact exchange contributions. 
However, in our constructions, MGGA_MS2 is a byprod- 
uct of MGGA_MS2h, which is fitted to the atomization 
energies of the G2/97 test set and the barrier heights 
of the BH42/03. Surprisingly, MGGA_MS2h performs 
reasonably well for the S22 set, much better than the 
other four hybrids, namely, PBEh, TPSSh, revTPSSh, 
and B3LYP. Note the improvement of MGGA_MS2h over 
the other four hybrids on the weak interactions is not a re- 
sult of sacrifice of accuracy for atomization energies and 
barrier heights as shown in Table [TTJ The good perfor- 
mance of MGGA_MS2h on the S22 set is then transfered 
to MGGA.MS2. 

When considering both Tables [XT] and |III| for the 

properties of atomization energies, barrier heights, and 
weak interactions, our three-parameter global hybrid 
MGGA_MS2h overall is the second best hybrid af- 
ter M06 that is heavily parameterized. Compared 
to MGGA_MS2h, our semilocal functional MGGA.MSl 



performs surprisingly well for atomization energies of 
molecules and barrier heights of chemical reactions al- 
though it contains just one empirical parameter. How- 
ever, the gain is obtained at the price of descriptions 
for weak interactions as demonstrated by the compari- 
son with MGGA_MS0 in Table HTH This shows the limit 
of tuning only the s~ dependence of the enhancement 
factor. Fortunately, by tuning the a— dependence addi- 
tionally, MGGA_MS2 improves over MGGAJVISO for the 
atomization energies significantly while remaining com- 
parable for the barrier heights and weak interactions. 



V. CONCLUSIONS 

By taking advantage of the understanding on the 
kinetic-energy-density dependence of MGGAs [17], we 
construct a global hybrid meta-generalized gradient ap- 
proximation (hybrid meta-GGA) with three empirical 
parameters (MGGA_MS2h), which has robust perfor- 



mances on the molecular systems for the properties of 
heats of formation, barrier heights, and noncovalent in- 
teractions. The derived underlying MGGAJVIS2 im- 
proves over the original MGGA_MS0 [T7] for heats of 
formation significantly while retaining good performance 
for barrier heights and weak interactions. By relaxing the 
parameter k of MGGA_MS0, we also obtained the one- 
parameter functional MGGA_MS1 that performs surpris- 
ingly well for heats of formations and barrier heights. 
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